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ABSTRACT 

We present axisymmetric numerical simulations of the solar interior, including 
the convection zone and an extended radiative interior. We find that differential 
rotation in the convection zone induces a toroidal field from an initially purely 
poloidal field. This toroidal field becomes unstable to the axisymmetric Tayler in- 
stability and undergoes equator-ward propagating toroidal field reversals. These 
reversals occur in the absence of a dynamo and without accompanying poloidal 
field reversals. The nature and presence of such reversals depends sensitively on 
the initial poloidal field strength imposed, with North-South symmetric reversals 
only seen at a particular initial field strength. Coupled with a dynamo mech- 
anism which regenerates the poloidal field this could be one ingredient in the 
sunspot cycle. 

Subject headings: solar physics, magnetohydrodynamics 
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1. Introduction 

One of the fundamental questions in solar physics is how the Sun generates its magnetic 
field. The most likely scenario is that the large-scale solar magnetic field is maintained 
by a hydromagnetic dynamo. This is supported by the remarkable regularity of the solar 
magnetic cycle. Over the course of the sunspot cycle bipolar pairs of magnetic regions 
appear at mid-latitudes over time, with new pairs appearing at lower and lower latitudes. 
These pairs have a distinct tilt in latitude (Joy's Law), with the leading sunspot (the one 
of the pair closest to the equator) having one polarity and the trailing sunspot having the 
opposite polarity. In the other hemisphere the polarities are reversed (Hale's law), so that 
the leading sunspots in the Northern hemisphere have the opposite sign of the leading 
sunspots in the Southern hemisphere. These sunspot pairs are associated with toroidal 
magentic field erupting to the surface, explaining their bi-polar nature. The leading sunspot 
polarity has the same polarity as the large scale poloidal field in that hemisphere and, after 
approximately 11 years, the sunspot cycle repeats, with the large scale magnetic polarity 
as well as that of the leading sunspot polarity reversing in each hemisphere reversed. The 
sunspot cycle itself is then said to have a period of approximately 11 years, while the global 
magnetic cycle has a period of approximately 22 years. 

Exactly how this dynamo process works, and produces these detailed observations 
of the sunspot and magnetic cycle remains the subject of intense research. Dynamo 
theorists usually envision a process whereby a large scale poloidal field is acted upon by 
differential rotation (most likely in the region of strongest shear, the tachocline), generating 
a strong toroidal field which then becomes magnetically buoyant and rises through the 
convection zone to produce sunspot pairs at the solar surface. Most research done on the 



solar dynamo has concentrated on solving mean field equations (Steenbeck et al. 1966 



Moffatt 1972), in which the induction of field due to correlations of fluctuating velocity 
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and magnetic field are prescribed by the so-called a effect and the differential rotation 
prescribed by the Q effect. These reduced models, with input from observations, can 
reproduce some of the features of the solar magnetic and sunspot cycle in the form of 



reversing toroidal fields which propagate toward the equator (Charbonneau 2005 Dikpati 



& Gilman 2008), albeit with tunable parameters. However, it has been more difficult to 



produce equatorward propagation of toroidal field in more realistic simulations which solve 



the full magnetohydrodynamics (MHD) equations (Glatzmaier 1985b[a Brun et al. 2004 



Brown et al. 2010). Early simulations (Glatzmaier 1985b) could produce a dynamo with 



oscillating toroidal field, but the period was usually much shorter than the 11 year sunspot 
cycle and toroidal field tended to propagate toward the pole. More recent simulations 



(Brown et al. 2010) have also been able to robustly produce dynamos, some of which show 



reversing toroidal field polarity, but again, the period is not solar like and are generally 



poleward propagating or stationary (Brown et al. 2010 Ghizaru et al. 2010). So it appears 



from these simulations that dynamos are robust, but dynamos which produce solar-like 
properties, such as equator-ward propagation of toroidal field are substantially more difficult 
to produce. 

While a hydromagnetic dynamo is the most likely explanation of the solar cycle, a 
magnetic oscillator is not completely ruled out. The dissipation time of a primordial large 
scale magnetic field in the Sun is ~ 10 10 years, therefore a decaying field superposed by an 
oscillating flow could possibly reproduce some of the solar observations. 

Most of the previous numerical simulations of dynamo action in the Sun and stars have 
focused on the convection zone, where it is thought that the rise and twist of magnetic field 



could produce a robust a effect. Spruit (1999, 2002) has argued that a dynamo could be 



driven in stellar radiative interiors requiring only differential rotation. In his picture the 



a effect could be generated by the Tayler instability. Braithwaite (2006) followed up this 
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analytic work with numerical simulations showing that a dynamo could be generated in 
stellar interiors, requiring only some small differential rotation. However, |Brun &: Zahn 



(2006) conducted similar numerical simulations of the solar radiative interior and found the 
instability but no regeneration of the mean poloidal field and therefore, concluded that no 
dynamo was present. Regardless of the presence of a dynamo, the presence of the Tayler 
instability in radiative interiors appears to be robust. 

The simulations presented here are axisymmetric and therefore, there is no dynamo. 
However, the simulations do reproduce some interesting features; namely, equatorward 
propagation of reversing toroidal field, in the absence of a dynamo or a reversing poloidal 
field. The rest of the paper is organized as follows. In section 2 we discuss the numerical 
model, in section 3 we discuss the evolution of the toroidal field, the role of the Tayler 
instability and the effect of the initial field strength. In section 4 we conclude with a 
discussion of the general applicability of this model to stellar interiors given its admittedly 
crude set up. 



Numerical Model 



We solve the full magnetohydrodynamic (MHD) equations in axisymmetric, spherical 
coordinates in the anelastic approximation. The numerical scheme is a combination of the 



equations outlined in Rogers & Glatzmaier (2005a) with the numerical scheme of Glatzmaier 



(1984). The most substantial differences between this model and the Glatzmaier (1984) 



model are that here we use a finite difference scheme in radius, our model is axisymmetric, 
and we solve the temperature equation rather than an entropy equation. For more details 
on the numerical method, see Appendix A. 
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2.1. Model Setup 

We solve the MHD equations in axisymmetric spherical coordinates in the radial range 
from 0.15 R to 0.95 R . The radiative region extends from 0.15 R Q to 0.71 R and the 
convection zone extends from 0.71 R to 0.90 R and we impose an additional stably 
stratified region from 0.90 R to 0.95 R . The reference state thermodynamic variables 



are taken from a polynomial fit to the standard solar model from Christensen-Dalsgaard 



et al. (1996). Therefore, the density ranges from approximately 60 gm/cm 3 at the bottom 
of the radiation zone to 10~ 2 gm/cm 3 at the top of the domain. The stable stratification is 
that given by the standard solar model and the super-adiabaticity of the convection zone 
is set to 10~ 7 . The grid resolution is 512 latitudinal grid points by 1500 radial grid points, 
with 500 radial zones dedicated to the radiative interior and 1000 zones dedicated to the 
tachocline and convection zone. 

Because this model is not fully three-dimensional (3D) we can not reproduce the 



Reynolds stresses required to set up the observed differential rotation (Thompson et al 



2003) in the convection zone and tachocline. We therefore, impose this differential rotation 



through a forcing term on the azimuthal component of the momentum equation: 
d , , (V • uUJt ((V x B) x B) t (2u x Q)^ 



dt r sin 9 r sin 9 pfir sin 9 r sin 9 



r sin 9 r 3 sin 9 t \ r sin 9 

where the last term on the rhs represents the forcing term. Q'(r,9) is the prescribed 
(observed) differential rotation relative to the constant Q and r is the e-folding time of the 
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forcing which we set to 10 4 Q The rotation profile is given by 

Sl'(r, 6) = to for r < 0.67^ 

n\r,6) = {A + B cos 2 6 + C cos 4 - 0)erf ( 2(r o " o ° 5 6 ^ 0) ) for 0.67 < r < 0.72 (2) 
fi'(r,0) = (A + Ecos^ + Ccos 4 ^-^) for r > 0.72 

where Q is the constant value set to 441nHz, A is 456nHz, B is -42nHz and C is -72 nHz 



(Thompson et al. 2003) 



Initially, we run a purely hydrodynamic model until convection is sufficiently established 
at which point we imposed a dipolar field with the form: 



so that initially all field lines close within the radiative interior, but overlap the bottom of 
the tachocline. We ran five models, which we label Models A,B,C,D, and E with initial field 
strengths, B Q , of 10G, 40G, 80G, 200G and 4000G, respectively, with all other parameters 
the same. Figure [T] shows a schematic of the initial setup and model, while Figure [2] 
shows the time-averaged azimuthal velocity as forced in Equation (2), along with the 
time-averaged meridional flow that results from this differential rotation profile (b). 



x We ran models with two other values of r (10 3 and 10 5 ) for a fraction of the total time 
of the models ran here and saw little difference in the meridional circulation or azimuthal 
flow variation. 
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3. Results 

3.1. Toroidal Field Evolution 

In the following we discuss the primary results from Model B (B a = 40G) and leave 
the discussion of the variation in field strength for Section 3.3. Because of the slight overlap 
of the imposed dipolar poloidal field with the tachocline and because of magnetic diffusion, 
it does not take long for the dipolar poloidal field to be stretched into toroidal field by the 
differential rotation in the tachocline and convection zone. Initially, this induces a toroidal 
field which has two signs each in the Northern and Southern hemisphere, see Figure [3} 
This structure can be understood by looking at the azimuthal component of the magnetic 
induction equation: 



dB$ rB d ^ d B^ sin 8 B d sin 8u e d B$ 

dt r dr r r dr r r 88 sin 8 r 88 sin 8 

+ B^u r h p + r/V 2 ^ - n ®\ a (4) 

r l sin 6 

The structure of the differential rotation gives du^/dr negative at high latitudes and 
positive at low latitudes and the first term on the rhs of (4) induces oppositely signed 
toroidal field at high and low latitudes, with asymmetry about the equator. This is seen in 
Figure |3^i. 

However, the first term on the rhs of (4), induction due to radial gradient in the 
rotation rate, is dominant for only a short time. The latitudinal toroidal field gradient 
becomes unstable (see Figure [6j section 3.2) and the dominant induction term becomes the 
fourth term on the rhs of (4), advection of toroidal field by latitudinal flow. This leads to a 
single dominant signed toroidal field in each hemisphere, asymmetric about the equator. In 
the following I will refer to positive toroidal field (white in figure 3) as the dominant sign 
in the Northern hemisphere, and negative toroidal field (black in figure 3) as the minority 
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field sign in the Northern hemisphere and conversely, in the Southern hemisphere. These 
particular signs arise naturally from the equator-ward directed ug coupled with a toroidal 
field which is stronger at the poles than the equator. Subsequently, and for most of the 
simulation, the advection of toroidal field by meridional flow (fourth term on the rhs of 
equation 4) is dominant within the tachocline (figure 3b). As can be seen from Figure |3]d to 
Figure ^jp a "reversal" takes place in which the dominant field is replaced by the minority 
field in the bulk of the convection zone and tachocline. This reversal proceeds by the 
advection of minority signed field toward the equator by meridional circulation at the base 
of the convection zone. 

The evolution in time of the reversal process is better demonstrated in a time-lattitude 
plot, as seen in Figure [4^,. Shown there is the toroidal field amplitude as a function of time 
and latitude. One can see the initial phase of two signs in each hemisphere giving way to 
a dominant sign in each hemisphere. Reversals subsequent to that shown in Figure [3] are 
also seen here. During the first reversal (around 7 x 10 7 s, or around 30 rotation periods) 
the northern and southern hemispheres are relatively in phase, however later reversals are 
stronger in the southern hemisphere and the coherence between north and south is not 
strictly maintained. As is seen in Figure |4| these reversals are seen both in the tachocline 
and at the top of the convection zone. We note that no reversal occurs in the radiative 
interior, nor does any reversal occur in the poloidal field. 

Except for the very initial phases the dominant term in the tachocline responsible 
for toroidal field evolution is the advection of toroidal field by the meridional circulation. 
The time averaged equator-ward latitudinal velocity at the base of the convection zone is 
3 x 10 3 cm/s, which would imply an advection time from pole to equator of ~ 3 x 10 7 s 
30 rotation periods). Looking at Figure [4] one can determine an approximate time for 
toroidal field to propagate from the pole to the equator to be ~ 1.5 x 10 7 — 2 x 10 7 s during 
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a reversal. This slight discrepancy is likely due to the difference between the time-averaged 
latitudinal velocity compared to the typical latitudinal velocity of the convective cells, 
which at the base of the convection zone is ~ 10 4 — 10 5 cm/s. This indicates that the bulk 
of the latitudinal transport is due to time-varying convection rather than the time-averaged 
meridional circulation]^] In animation of the images shown in Figure |3^i-d one can see that 
the toroidal field is mostly advected by the convective eddies, therefore leading to a slightly 
more rapid equator-ward propagation than would be implied by the mean latitudinal flow. 

Over the duration of the simulation, the polodial field does little else than diffuse and 
be advected by the overlying convection zone, as can be seen in figure 3e-h. As this model 
is axisymmetric there can clearly be no regeneration of the poloidal field and hence, no 
dynamo. Therefore, one would expect the reversals of field sign to be due to fluctuations of 
the azimuthal velocity superimposed on the imposed differential rotation. In Figure [5] we 
show the azimuthal velocity as a function of time and latitude. One can see the general 
behavior of the (imposed) observed differential rotation with fast equator (prograde flow) 
and slow poles (retrograde flow), this is maintained throughout the simulation. However, 
one can also see weak signs of prograde flow at both poles, intermittent in the tachocline 
but persistent in the convection zone. Likewise one can see bands of prograde azimuthal 
flows starting at mid-latitudes with slight propagation toward the equator, although on 
much shorter timescales than the magnetic reversals, contrary to observations. This short 
timescale structure with equator- ward propagation is also seen in the Figures 6-8, in the 

2 Note that the convective velocities in this simulation are higher than that expected in 
the solar convection zone. This is because our enhanced thermal diffusion leads to a larger 
heat flux through the system, requiring a larger convective velocity to carry it, given the 
same temperature gradient. In these simulations these effect leads to an enhancement of the 
flux by 10 6 and hence, convective velocities enhanced by ~ 100. 
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Tayler instability panel and, at some level, in the ratio of toroidal to poloidal field indicating 
that there is some interplay between magnetism and azimuthal flow variations, we are 
currently pursuing this interaction. 

Because of the lack of obvious variation in the azimuthal flow on the same timescales 
of the magnetic field it appears that the reversals must be due to a local change in sign 
and amplification in toroidal field followed by advection. This can happen if the persistent 
minority sign seen at high latitudes is amplified and becomes dominant in the regionF] 



3.2. Tayler Instability and toroidal field evolution 



There has been significant discussion about the possible role of magnetic instabilities 



in stellar interiors (Wright 1973 Tayler 1975 Acheson 1978 Spruit 1999 Parfrey & Menou 



2007). Although it is generally believed that convective or turbulent motions provide 



the required a effect to regenerate poloidal field from toroidal field, it was suggested by 



Spruit (2002), that convection is not needed because this effect could be provided by an 



instability in the field itself. He argued that the most relevant instability in stellar interiors 



was likely to be the Tayler instability. Parfrey & Menou (2007) argued that the diffusive 



magneto-rotational instability (MRI) was likely at high latitudes in the tachocline and that 
while this region overlaps with the region of Tayler instability the MRI is likely to be more 
dynamically important because of its exponential growth. They argued that this instability 
could disrupt magnetic fields, thus explaining the lack of sunspots at high latitude. 

Here we find that the toroidal field reversals observed in Figure [4] are consistent with 
times when the toroidal field is unstable to the axisymmetric Tayler instability as described 



3 This may only be possible with prograde poles, which we see here but are not necessarily 
realistic. 
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by the simple instability criterion provided in Goossens et al. (1981) and Spruit (1999), 
namely: 



d / , B 



2 



cos "m{ ln ^ ) ) >0 (5 > 

In Figure [6Jd we show the lhs of the inequality of (5). One can clearly see that 
reversals of the field are coincident with a Tayler instability that proceeds even at low 
latitudes. The instability does not appear in the bulk of the radiative interior where the 
toroidal field strength is substantially lower due to weaker differential rotation and the 
(stabilizing) poloidal field is correspondingly stronger. The instability is strongest in the 
tachocline and is present, but much less coherent in the convection zone. The bulk of 



the time the instability is confined to higher latitudes as expected (Parfrey & Menou 



2007). In Figure [6ji, which shows the ratio of the local toroidal to poloidal field energy 
at high latitudes, one can see an almost oscillatory behavior of the instability. At these 
high latitudes the dominant toroidal field grows due to differential rotation, once the local 
ratio reaches a value around 10-20, the instability disrupts the field growth and magnetic 
flux destruction ensues. Normally, the field decay is weak and growth due to differential 
rotation is re-established shortly thereafter. However, occasionally, field destruction caused 
by the instability is more severe and the toroidal field energy drops precipitously. During 
these times the minority signed field can become dominant and be advected equator-ward. 
Therefore, times of field reversals as seen in Figure [6j [7] and [8] are associated with minima 
of the toroidal field energy relative to the poloidal energy and these minima are induced by 
the Tayler instability. 

Indeed, one can see in Figure [6]i that the slight North-South asymmetry in reversal 
coincides with a slight discrepancy in timing between flux destruction in the Northern and 
Southern hemispheres. This can be seen in all of the composite figures (6-8), when flux 
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destruction occurs there is a reversal, when there is weak or no destruction, no reversal 
ensues. 

The stability condition expressed in (5) is for a purely toroidal field. It has been 



shown (Wright 1973; Braithwaite & Spruit 2004) that mixed field configurations are more 



stable. Therefore, in addition to the condition expressed in (5) the toroidal field must also 
overcome the stabilizing effect of the poloidal field. This is what we see in the bottom 
panels of Figures [6j [7] and |8j Substantial flux destruction only occurs when the criterion 
in (5) is met, along with the ratio of local toroidal to poloidal magnetic energy reaching 
some critical value. The critical value appears to be around 10-20, but clearly is not strict 
and other factors, such as local turbulent dissipation could become relevant. It should be 
noted that the non-axisymmetric Tayler instability and MRI are likely more relevant in real 
stars, this could lead to more regular flux destruction and possibly more regular reversals, 
although whether they would proceed like the ones here is unclear. 



3.3. Dependence on Initial Field Strength 

In the preceding section we concentrated on Model B and found that it showed 
equatorward propagating toroidal field reversals. These reversals were precipitated 
by a Tayler instability and were sometimes coherent between Northern and Southern 
hemispheres, but not always. We ran four additional models with varying initial field 
strengths. All of the models, except Model A, show some evidence of reversal, however, 
only Model B shows reversals which are coherent in Northern and Southern hemispheres. 
Furthermore, only Model B and E show reversals of any substantial duration. Model A, 
which shows no reversal does show dominant toroidal field which weakens and then increases 
in time, but this is not accompanied by amplification and advection of minority field. 
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Figures 7 and 8 shows the toroidal magnetic field, Tayler instability criterion and 
ratio of toroidal to poloidal magnetic field energy for Model E and C, respectively. One 
can see the initial two signed field structure giving way to a single dominant field in each 
hemisphere, similar to Model B. Other equatorial propagting field reversals are also seen, 
particularly in the Southern hemisphere for Model E and the Northern hemisphere for 
Model C. One should note that the local ratio of the toroidal to poloidal field strengths 
in both of these models is similar to that seen in Model B, despite the vast difference in 
initial poloidal field strength; again implying the ratio of the field strengths plays some 
role in determining the growth of the Tayler instability. It is unclear why Model B shows 
symmetric reversals, while the others do not, it is also not clear why the various reversals 
are so different both in latitudinal excursion and in duration. These simulations indicate 
that weakening of the dominant toroidal field by the Tayler instability is robust, however, 
accompanying solar-like reversals may be possible only in limited parameter space of 
magnetic field strength. The dynamics that dictate that parameter space are presently 
unclear. 



4. Discussion 

We have shown that a poloidal field, initially confined to the radiative interior and 
acted on by differential rotation of the kind observed in the Sun can lead to reversals and 
equator-ward propagation of toroidal field. The process is initiated by differential rotation 
acting on an initially poloidal field to induce toroidal field. When both the basic Tayler 
instability criterion expressed in (5) and the local ratio of toroidal to poloidal field energy 
reaches some critical value, instability and rapid flux destruction ensues. During this time 
the minority field sign can become dominant and be transported toward the equator by 
meridional circulation. Once the minority field decays and is advected away the dominant 
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field remains and the process repeats. 

The relevance of this process in the Sun or any star are unknown. First, the toroidal 
field reversals do not appear periodic, nor completely symmetric about the equator, although 
this depends on the initial field strength. Second, the non-axisymmetric Tayler instability 
or the weak field MRI is likely dominant in the solar radiative interior and tachocline. 
However, any instability that could lead to the amplification of the minority signed toroidal 
field could feasibly mimic what we see here. In fact, a more robust instability may lead to 
more regular reversals. Third, the high latitude fluctuations of azimuthal velocity which 
initiate this process, while out of sight of helioseismology, may be unrealistic. Finally, 
this process clearly indicates a high-latitude branch of toroidal field which is equator-ward 
propagating, contrary to the observations. 

While these are all serious issues which complicate the applicability of this model to 
the Sun, we think there are substantial relevant overlaps with the Sun that make this 
model interesting. First, this is the first model which solves the full MHD equations 
with convection and rotation and gets equator symmetric equator-ward propagating and 
reversing toroidal fielcQ and, interestingly, it does so without the operation of a dynamo and 
without corresponding poloidal field reversals.. We also find that these reversals depend 
on the initially imposed field strength, and it is possible that there is a preferred field 
strength which more adequately represents observations. Of course, we would like to run 
more models to confirm this, but computational resources are limited. Second, we have 
shown that the axisymmetric Tayler instability, likely to be operating in stellar interiors, 



^ote, however, that 



Mitra et al. 



(????) have simulated equator- ward propagating re- 
versing fields using forced MHD turbulence with imposed helicity. Those models, implicitly 
including rotation and convection through helicity have produced much more regular rever- 
sals. 
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can act to destroy the dominant toroidal field, allowing the weaker minority field to be 
advected equator-ward by the meridional circulation at the base of the convection zone. 
This instability seems to act when the basic criterion expressed in (5) is met combined 
with when the local ratio of toroidal to poloidal field energies reaches some critical value. 
Third, we have shown that meridional circulation can effectively act as a conveyor belt in 
advecting, at least in our 2D geometry, minority field toward the equator at the base of 



the convection zone; a requirement for many flux-transport dynamo models (Charbonneau 



2005 Dikpati et al. 2007). 
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Appendix A: Numerical Method 

In this model we solve the full set of axisymmetric MHD equations in the anelastic 
approximation. 



V • (pu) = 0. 
V-B = 



(1) 
(2) 



-^ + (u-V)u = -VP-C^f + 2(ux ft) + -(J xB) 



+ v v^u + -v(v-u; 



(3) 



dT 
~dt 



+ (u- V)T 



(dT 
Tr - (7 - l)Th p 



+ (7 - l)Th p u r + 7 k[V 2 T + (h p + K 



dB 

~dt 



V x (u x B) + r/V 2 B 



dT 
dr 



(4) 



(5) 



Equations 1 and 2 ensure the mass and magnetic flux are conserved. In Equation 3, the 
momentum equation, g is gravity, Q is the rotation rate and C denotes the co-density 



(Braginsky & Roberts 1995 Rogers & Glatzmaier 2005a), defined as: 



T \ gp or 



(6) 
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where T and p are the reference state temperature and density, and dT/dr is the reference 
state temperature gradient. T, p and P are the perturbation temperature, density and 
pressure, while u is the velocity with components u r ,ue and u^. 

Equation 4 represents the energy equation written as a temperature equation, where 7 
is the adiabatic index (which we specify to be 5/3) and h p and h K are the inverse density and 



thermal diffusivity scale heights, respectively. As in Rogers & Glatzmaier (2005a), we utilize 



the temperature equation in this formulation of the anelastic equations, as opposed to the 



entropy formulation commonly used (Glatzmaier 1985b). This allows us to better specify 



the super- and subadiabaticity of regions, which can be seen in Equation 4, where the first 
term on the right hand side (rhs) represents the difference between the reference state 
temperature gradient and the adiabatic temperature gradient and allows us to represent 
strongly subadiabatic regions in addition to convection zones. The last term involving Q 
represents the physics included in the standard solar model, which is not included in this 
model, that maintains the initial reference state temperature gradient. The sum of the 
last two terms in Equation 4 account for the total reference state flux through the system, 
which is set to zero, so that the reference state temperature is time-independent. 

Equation 5 is the magnetic induction equation, in which n represents the magnetic 
diffusivity. The thermal diffusivity k is specified to have the solar profile but multiplied by a 



constant for numerical reasons, as in Rogers & Glatzmaier (2005a). The viscous diffusivity 



is inversely proportional to the density and the magnetic diffusivity is constant: 

3 

K - K 16aT 



'■3p 2 kc p 



V = (7) 
P 



T] = const 

where a is the Stefan-Boltzman constant, k is the opacity and c p is the specific heat at 
constant pressure. In the simulations presented here K m is 10 6 , v m is 5 x 10 11 and rj is 
5 x 10 12 . Typical values for various parameters are shown in Table 1. 
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Parameter 


Symbol 


Sun 


Simulation 


thermal diffusivity 


K 


10 7 


1.1 x 10 13 


magnetic diffusivity 


V 


10 3 


5 x 10 12 


viscous diffusivity 


V 


30 


2.5 x 10 12 


Rayleigh number 


gD 4 AVT sup /T bcz un 


10 23 


10 6 


Prandtl number 


V j K 


3 x 10~ 6 


0.22 


magnetic Prandtl nb 


v/j) 


3 x 10~ 2 


0.5 


rotation frequency 


Q 


2.7 x 10~ 6 


2.7 x 10~ 6 


Froude number 


(n/N) 2 


2 x 10" 6 


2 x 10~ 6 


Ekman number 


v/2VLR 2 


2 x 10~ 15 


2 x 10" 4 



Table 1: Values for various parameters in the simulations. The values listed are those at 
the base of the convection zone/top of the radiation zone. Diffusion times are calculated 
using the depth of the radiation zone and the diffusion values at the base of the convection 



zone/top of the radiation zone. The Rayleigh number is calculated as in Rogers & Glatzmaier 



(2005b) and here we calculate the solar value assuming a super-adiabaticity of 10~ 7 . Models 
were run approximately 2 x 10 8 s and each required approximately a few hundred thousand 
processor hours. 
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The method here is similar to that outlined in Glatzmaier (1984) in which the mass 



and magnetic flux are written in terms of poloidal and toroidal components: 

pu = V x V x (Wr) + V x (Zf) 
B = V x V x (Br) + V x (Jf) 



(8) 
(9) 



The independent variables W, Z, B, J and the thermodynamic variables temperature, T, 
and pressure P, are then expanded in spherical harmonics so that: 



pu r 



pug 



l -Y,Kl + mi{r,t)Pi 
i 

in ft ^— — ' f)r f)ft 



r sin 9 



^£-^(r,t)sin0 



r sin 9 



09 

m 

09 



and similarly for the magnetic field 

B r = \j2^ l+1 ^ Bl ^^ Pl 



1 ^ dBijryt) . 0P t 

r sm 9 z — ' Or 09 

i 



— ?— Y)-J,(r,*) sin# 
r sm 9 z — ' 



op 

09 



(10) 



in 



;i2) 



(13) 

(14) 
(15) 



where 9 is the co-latitude, I is the spherical harmonic degree and Pi are the Legendre 
polynomials. Temperature and pressure are similarly expanded so that: 



T = Y,T l (r,t)P l 
i 

p = Y, v M p i 



(16) 
(17) 



-23- 



The radial functions Wi(r,t), Zi(r,t), B[(r,t), Ji(r,t), Pi(r,t) and Ti(r,t) are discretized 
on a non-uniform, second-order finite difference grid in radius. This is done to allow us 
the flexibility of enhanced resolution around the tachocline and overshoot layer, without 
having the burden of insufficient resolution at the tachocline, required by the Chebyshev 
grid or alternatively, enhanced resolution at the bottom of the radiation zone (as would be 
necessary with stacked Chebyshev grids). 

We then solve the radial component of the momentum equation for Wi, the radial 
component of the curl of the momentum equation for Zi and the full divergence of the 
momentum equation for Pi. The equation for Pi is then a diagnostic equation, unlike the 



Glatzmaier code Glatzmaier (1984) which solve the horizontal divergence of the momentum 
equation. This is done to avoid third order derivatives, which would complicate boundary 
conditions when using the finite difference scheme employed here. We impose the stress-free 
and impermeable boundary conditions on the velocity. 

For the magnetic terms we solve the radial component of the magnetic induction 
equation for Bi and the radial component of the curl of the magnetic induction equation 
for Ji. The temperature and pressure equations are solved directly for 7] and Pi. The 
boundary conditions on temperature are constant temperature perturbation, on velocity we 
force stress-free and impermeable conditions, the magnetic field is matched to an internal 
potential field at the bottom of the domain and an external potential field at the top of the 
domain. The current density, Ji is zero on the top and bottom boundary. 

Timestepping is accomplished by an explicit Adams-Bashforth method for the 
non-linear terms and an implicit Crank-Nicolson method for the linear terms. The code 
is parallelized using Message-Passing- Interface (MPI). Each model was run approximately 
1 — 2 x 10 8 s and took approximately a few hundred thousand processor hours. 



Fig. 1. — Model schematic. Radiative interior extends from 0.15R Q to 0.71 R and con- 
vection zone lying from O.71R to O.9OR . After convection is established a dipolar field is 
imposed in the radiative interior. Waves of mixed gravity-Alfven type exist in the radiative 
interior. 
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Fig. 2. — Time averaged differential rotation for Model B is shown in (a) with red repre- 
senting flows faster than the mean and blue representing flows slower than the mean. Red 
overlaid lines represent the extent of the initially imposed tachocline. The figure shows the 
forced differential rotation as described in the text, plus any fluctuations averaged in time. 
The amplitudes in the color bar are in Hz. Time averaged meridional circulation for Model 
B (b), white colors represent flow towards the south pole, while dark colors represent flow 
towards the north pole. One can then clearly see poleward flow at the top of the convection 
zone, and equator-ward flow at the base of the convection zone. Amplitudes in the color 
bar are quoted in cm/s. White overlaid lines represent the extent of the initially imposed 
tachocline. Both time-averages are over the entire simulation (2 x 10 8 s). 
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(o) 8 x 10 6 s (b) 7 x 10' s (c] 8 x 10' s (d) 1 8 s 




(e) (f) (g) (h) 




Fig. 3. — Time snapshots of the toroidal field (a-d) and poloidal field (e-h). In (a-d) white 
represents positive toroidal field, while black represents negative. Initially, the radial gradient 
of the differential rotation dominates the induction of toroidal field, this as seen as two signs 
of toroidal field in each hemisphere (a). However, quickly the advection of toroidal field by 
the meridional flow becomes the dominant term controlling toroidal field evolution and one 
sign in each hemisphere remains. A toroidal field reversal in the tachocline and convection 
zone is seen between (b) and (c), see corresponding time in Figure |4j Poloidal field is seen in 
(e-h) and is diffusing outward and decaying in time. The colormaps in the various images are 
slightly different because of the steady increase in toroidal field strength and steady decrease 
in poloidal field strength. Peak toroidal field in this model is approximately ±50G. Both 
poloidal and toroidal field effectively connect to the convection zone despite the meridional 
circulation induced by the differential rotation. 
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Fig. 4. — Toroidal magnetic field as a function of time and latitude, in the tachocline 
(bottom) and at the top of the convection zone (top). Toroidal field reversals with equator- 
ward propagation are seen. The Northern hemisphere shows a weaker second reversal and 
no third reversal. 



-28- 




Fig. 5. — Azimuthal velocity as a function of time and latitude, in the tachocline (bottom) 
and in the convection zone (top). One can see some signs of short period torsional oscillations 
which propagate toward the equator in time. These are symmetric about the equator and 
relatively periodic. However, their timescale is more like that of convective eddies rather 
than magnetic field reversals, unlike the Sun. 
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Fig. 6. — Toroidal magnetic field as a function of time and latitude within the tachocline, 
with color scaled at ±30G (a), Tayler instability criterion as described in Equation 5 with 
color scaled at ±1 (b), ratio of toroidal to poloidal magnetic field energies, with color scaled 
0-300 and white representing high amplitude ratios (c) and ratio of toroidal to poloidal 
magnetic field energies at high latitudes (20°) in the Northern hemisphere (solid line) and 
Southern hemisphere (dotted line) (d). The right hand panels show the same values just 
zoomed over the first reversal Deriod. 
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Fig. 7. — Same as Figure [6j but here showing Model E (B = AkG) and neglecting the 
zoomed panels. Because this model has a larger initial field strength, (a) is scaled at ±50G, 
others are scaled the same as Figure |6j 
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Fig. 8. — Same as Figure [6j but here showing Model C (B = 80(2) • Because this model has 
a weaker initial field strength, (a) is scaled at ±3. 



